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ABSTRACT 

We study the heating of charged test particles in three-dimensional numerical simulations of weakly 
compressible magnetohydrodynamic (MHD) turbulence ("Alfvenic turbulence"); these results are rel- 
evant to particle heating and acceleration in the solar wind, solar flares, accretion disks onto black 
holes, and other astrophysics and heliospheric environments. The physics of particle heating depends 
on whether the gyrofrequency of a particle ilg is comparable to the frequency of a turbulent fluctuation 
Lo that is resolved on the computational domain. Particles with ^l^ ^ uj undergo strong perpendic- 
ular heating (relative to the local magnetic field) and pitch angle scattering. By contrast, particles 
with Lj undergo strong parallel heating. Simulations with a finite resistivity produce additional 

parallel heating due to parallel electric fields in small-scale current sheets. Many of our results are 
consistent with linear theory predictions for the particle heating produced by the Alfven and slow 
magnetosonic waves that make up Alfvenic turbulence. However, in contrast to linear theory predic- 
tions, energy exchange is not dominated by discrete resonances between particles and waves; instead, 
the resonances are substantially "broadened." We discuss the implications of our results for solar and 
astrophysics problems, in particular the thermodynamics of the near-Earth solar wind. This requires 
an extrapolation of our results to higher numerical resolution, because the dynamic range that can 
be simulated is far less than the true dynamic range between the proton cyclotron frequency and the 
outer-scale frequency of MHD turbulence. We conclude that Alfvenic turbulence produces significant 
parallel heating via the interaction between particles and magnetic field compressions ("slow waves"). 
However, on scales above the proton Larmor radius Alfvenic turbulence does not produce significant 
perpendicular heating of protons or minor ions (this is consistent with linear theory, but inconsistent 
with previous claims from test particle simulations). Instead, the Alfven wave energy cascades to 
perpendicular scales below the proton Larmor radius, initiating a kinetic Alfven wave cascade. 
Subject headings: turbulence - MHD - plasma - particle acceleration - solar wind - accretion disks 



L INTRODUCTION 

The heating and acceleration of particles by magne- 
tohydrodynamic (MHD) turbulence is believed to play 
a critical role in phenome na as diverse as the solar 
wind and solar corona (e.g., Cranmer & van Ballegooi-i 
2005|), accretion disks in active galactic nuclei (e.g.. 



Cjuataert fe Gruzinov] 19991, and the confineme nt and 



re-accelera t ion of cosmic-rays in th e Galaxy (e.g., Chan 



dran 2000 Yan & Lazarian 2002|. This paper focuses 
on the physics of particle heating by weakly compressible 
MHD turbulence. Weakly compressible MHD turbulence 
consists of nonlinearly interacting Alfven and slow mag- 
netosonic waves; since the Alfven wa ves dominate the 
dynamics fColdreich & Sridhar 19951, we will often re- 
fer to such turbulence as Alfvenic turbulence. Measure- 
ments of magnetized t urbulence in laboratory plasmas 



and in the solar wind (Bale et al. 



2005 Sahraoui et al. 



2009|), as well as analytic models Shebalin et al. 1983 



Higdon 1984 



Goldreich fc S ridhar T^MT and numerical 



simulations (Cho & Vishniac 2000: Maron & GoldreichJ 



2001 1 demonstrate that most of the energy in Alfvenic 
turbulence cascades to small scales perpendicular to the 
magnetic field. This strongly infl uences how Alfvenic tur- 
bulence couples to particles (e.g., Quataert 1998 Leamon 
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etal]p99 l. 

i'he solar wind and solar corona provide a particularly 
rich source of data on both the properties of MHD tur- 
bulence (e.g., the magnetic and electric field power spec- 
tra) and on the possible thermodynamic outcome of par- 
ticle heating and acceleration by such turbulence (e.g., 
the proton, electron, and minor ion temperatures and/or 
distribution functions) . A detailed review of these obser- 
vational results and their theore tical interpret ation is not 
critical for this paper (see, e.g., Marsch|[2006 1 . However, 
one result that is particularly germane is the fact that 
the proton distribution function is typically anisotropic 
with respect to the local magnetic field: Tl — 0.9 T|| 
in the solar wind at ~ 1 AU, although t he sign of the 
anisqtropy depends on th e wind speed (Kasper et al. 



2002 Hellinger et al.|2006 1 . In the fast solar wind, minor 



ions such as O VI have T±_ ^ Ty, with Tj_ /r|| in creasing 
with the charge to mass ratio of the ions (Marsch 2006 ) . 
An outstanding problem is whether these measurements 
can be accounted for by heating by anisotropic Alfvenic 
turbulence - this is a particularly important question 
given the growing body of evidence that the turbulent 
fluctuations in the solar wind at '--^ 1 AU are consistent 



with anisotropic Alf venic turbulence (Bale et al. 2005 



Sahraoui et al]|2009 1 



i'he question of how Alfvenic turbulence heats and ac- 
celerates particles is equally pressing in other astrophys- 
ical contexts. For example, in a particular class of mod- 
els for accretion onto black holes. Coulomb collisions be- 
tween electrons and ions are too slow to maintain thermal 
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equilibri um, resulting in a two temperature plasma with 
Ti ^ Te dRees et al.||l982|. Weakly compressible MHD 



turbulence is present m accretion disks because of the 
nonli near saturation of th e magnetorotational instability 
( |Balbus fc Hawley 1998). In two-temperature accretion 
disk models, the radiation from the accreting plasma is 
produced largely by the (hghter) electrons, and is thus 
sensitive to the details of how the turbulent energy is 
thermalized as heat. 

There is an extensive body of work studying the heat- 
ing and acceleration of (test) particles by plasma waves 
using "quasilinear" theory. The key assumptions of this 
theory are that the orbits of the particles are given by un- 
perturbed helical motion around magnetic field lines and 
that the plasma waves are long-lived; as a result, energy 
exchange between waves and particles occurs only at dis- 
crete resonances (e.g.,|Stix 1992 and references therein; 
see To the extent that MHD turbulence can be ap- 
proximated as a superposition of roughly linear waves, 
the same resonances should describe how magnetized tur- 
bulence couples to the underlying plasma. This view- 
point dominates the literature on par ticle heating and 



acceleration by MHD t urbulence (e.g., [Miller fc Roberts 
T995l|Quataert. l 99 8: Leamon et al.|1999||Chandran|2nUnf 
Cranmer fc van Ballegooijen||2003y ! 

In this paper, we go beyond linear theory by directly 
simulating the orbits of charged test particles in three- 
dimensional numerical simulations of weakly compress- 
ible MHD turbulence. Our numerical simulations cap- 
ture the physics of wave-particle interactions (e.g., cy- 
clotron resonances) without making many of the simpli- 
fying assumptions of standard linear theory. In addition 
to wave-particle resonances that can in principle operate 
throughout the inertial range in a turbulent plasma, cur- 
rent sheets can form on small scales; these current sheets 
represent an additional mechanism for particle heating. 
In order to isolate the effects of current sheets, we calcu- 
late test particle heating in simulations with and without 
an explicit resistivity; wc note up-front, however, that 
our MHD simulations probably do not correctly capture 
how reconnection heats particles (see Sjs]) . Previous work 
on test particle hea ting in numer ical simulations of MHD 
turbulence (Dmitruk et al. 2004 1 concluded that the pref- 
erential perpendicular heating of ions in the solar wind 
summarized above could readily be accounted for. We 
disagree with this conclusion for reasons discussed in SjSj 

The test particle method used here has significant 
strengths, but also some weaknesses. Our grid-based 
MHD code can resolve a reasonably large inertial range 
in a computationally efficient manner with test parti- 
cles spanning many orders of magnitude in gyrofrequency 
and/or velocity, thus accurately capturing much of the 
relevant physics in a single calculation. The downside 
of the test particle approximation is that efficient par- 
ticle heating or acceleration would in reality modify the 
turbulent cascade, an effect that is not captured by our 
approach. Fully electromagnetic particle-in-cell (PIC) 
calculations of Alfvenic turbulence are not feasible with 
current computing power. An alternative approach, gy- 
rokinetics, self-consistently captures the physics of low 
frequency kinetic turbulence, and can be used to study 
the combined problem of particle heating and its effect 
on the turbulent cascade (Howes et al.^ ,2008b). Gyroki- 
netic simulations arc, however, computationally intensive 



and order out higher frequency dynamics such as the cy- 
clotron resonances. The test particle method used here 
is thus a complementary approach that can be used to 
rapidly and accurately explore the wave-particle interac- 
tions and reconnection physics that occur on large-scales 
in MHD turbulence. 

The remainder of this paper is organized as follows. ^ 
reviews the linear theory predictions for particle heating 
by Alfvenic turbulence; these provide a useful framework 
for interpreting our numerical results. S|3]outlines our nu- 
merical methods, with some of the details given in the 
Appen dix. We present the results of a fiducial simula- 
tion in |4.1| along with our interpretation of the resulting 
particle heating and its relation to linear theory; §^4.2 



4.5 explore the effects of varying the resistivity, magnetic 
field strength, and particle distribution function. Finally, 
in S|5] we summarize our results and discuss their impli- 
cations, focusing on the near-Earth solar wind. 

2. ANALYTIC EXPECTATIONS 

Weakly compressible Alfvenic turbulence can be 
viewed as a collection of nonline arly interacting Alfven 
and slow magnetosonic waves (Goldreich fc Sridhar 
19951. We focus largely, but not entirely, on plasmas 
with P > 1, so that the linear dispersion relation of both 
waves is given by w ~ |fc|| where uj is the frequency of 
the wave, fcy is its wave-vector along the local magnetic 
field, and va is the local Alfven speed. In quasilinear 
theory, these waves exchange ener gy with pa rticles only 
at discrete resonances, when (e.g., Stix||1992 1 

LO — — (1) 

where is the particle's velocity along the magnetic 
field, is the particle's cyclotron frequency, and n is an 
integer. When equation M is satisfied for a linear wave, 
there is phase coherence averaged over long timescales 
^ w"^, and thus energy can be exchanged between the 
wave and the particles. In strong Alfvenic turbulence, 
however, nonlinear interactions transfer energy from one 
scale to another on a timescale comparable to the lin- 



ear wave period (^ 



). It is thus unclear whether 



standard quasilinear theory will adequately describe the 
heating of particles by strong Alfvenic turbulence. At a 
minimum, we expect the discrete resonances in equation 
([T]) - which show up as delta functions in linear theory - 
to be significantly "broadened" due to the rapid decorre- 
lation of the ph ases of waves in strong turbulence (e.g., 
Chandran|[2000l ). 

T'he n = resonance in equation ([T]) generally applies 
to low frequency fiuctuations having ui <^ ^l. In this 
case resonance occurs when lo = i.e., when the 

parallel velocity of a particle is comparable to the parallel 
phase speed of the wave, so that particles "surf" along 
with the wave. In the presence of such low frequency 
electromagnetic fluctuations, the magnetic moment of a 
particle, defined here as^ 



M = 



B 



(2) 



is an adiabatic invariant and remains approximately con- 
stant {B = \B\ is the magnitude of the magnetic field 
strength and u± is the gyration/perpendicular velocity). 

* In equation |2|, we have dropped the conventional factor of 
m/2, since the masses of the particles we consider are arbitrary. 
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It is straightforward to show that in the hniit uj <^ fl, 
the equation of moti on for a charged particle simplifies 
to ( |Achterberg||1981[ ) 

du\\ a , , 

/iViiB (3) 



dt 



and 



du±^ 
"df 



2B 



dB 
~dt 



(4) 



where E\\ is the parallel electric field and V|| is the gra- 
dient along the local magnetic field; note that equation 
Q refers to the magnitude of the perpendicular velocity 
and so does not contain the cyclotron motion. 

Equation ([3| demonstrates that the ri = resonance 
contains two physically distinct heating mechanisms, 
that due to parallel electric fields (Landau damping) and 
that due to the fi\7B force (transit-time damping); for 
simplicity, however, we will refer to this low-frequency 
resonant interaction as the "Landau resonance." In 
MHD, linear Alfven waves have 5E^^,SB = and thus 
cannot interact with particles via the Landau resonance. 
Alfven waves do develop SEh and/or SB ^ on small 



scales where MHD breaks down (e.g., Quataert 



19981, 



but this is not captured by our turbulence simulations 
which utilize MHD. 

Unlike the Alfven wave, the slow wave has 5B 7^ in 
MHD ; analytic studies of particle heating by anisotropic 
Alfvenic turbulence have thus predicted that the slow 
wave should produc e significant hea ting via the Lan- 
dau resonance (e.g., 
this is predicted to 



Blackman 19991. In linear theory, 
3e entirely parallel heating, i.e., it 
should increase the parallel temperature, but not the 
perpendicular temperature, of the particles. This fol- 
lows directly from equation Q by considering a sin- 
gle Fourier component with 6B cx g-«'^*+«k r^ -^^ which 
case du±/dt oc [w — By construction, however, 

Lo = at the Landau resonance and so du±/dt = 0. 

li uj ^ il, i.e., if the fields vary significantly on the 
timescale of a particle's cyclotron motion, strong "cy- 
clotron resonance" may occur; this corresponds to n — 
±1 in equation Q.^ Cyclotron resonance occurs when 
the perpendicular electric force due to a wave remains in 
phase with the rotating cyclotron motion of a particle. 
Because a; f2, adiabatic invariance of /i is violated and 
cyclotron resonance leads to strong pitch angle scattering 
and perpendicular heating. 

The Alfven wave component of anisotropic MHD tur- 
bulence is a transverse linearly polarized wave that can 
undergo cyclotron resonance with both negatively and 
positively charged particles. By contrast, anisotropic 
slow waves are largely longitudinal and thus do not un- 
dergo strong cyclotron resonance. To estimate the con- 
ditions required for cyclotron resonance in anisotropic 
MHD turbulence, we use Goldreich & Sridhar's (1995) 
critical balance conjecture that relates the typical paral- 
lel (fc||) and perpendicular (fc_L) wavevectors of turbulent 



fluctuations: fcii 



i,2/3. 1/3 



where k^- is the outer 



scale of the turbulence. This scale-dependent anisotropy 
of incompressible MHD turbulen ce has been confirmed 
in a number of numerical studies ( Cho & Vishniac 2000} 
Maron fc GoldreichpOOl I. Together with the dispersion 



relation for Alfven waves in MHD, lo ~ |fc|||w^, this im- 
plies that cyclotron resonance occurs when 

^A-un^\=±n. (5) 



1,2/3,1/3 
"-i "-mm 



Equation ^ will be useful for interpreting some of the 
numerical results described later in this paper. 

To summarize the above discussion, quasilinear theory 
predicts that anisotropic Alfvenic turbulence resonantly 
couples to particles in two ways: (1) parallel /iV||i3 
heating of particles by Landau-resonant slow modes, (2) 
pitch-angle scattering and perpendicular heating of par- 
ticles by cyclotron-resonant Alfven waves. Quasilinear 
theory does not account for possible heating of particles 
at current sheets, but this will naturally be a feature of 
our resistive MHD simulations, in addition to the wave- 
particle resonances reviewed here. 

3. NUMERICAL METHODS 

Our simulations focus on coUisionless test particles 
evolving in isothermal subsonic MHD turbulence. Our 
approach involves two distinct computational phases. 
First, the macroscopic fields are evolved according to the 
MHD equations. Then the test particles' positions and 
velocities are updated, using the macroscopic fields to 
compute the appropriate Lorentz force. 

3.1. The MHD Integrator 
To compu te the macroscopic fields , we use the Athen a 



MHD code ( [Gardiner fc Stone 
Calculations are done on a uni 



2005 Stone et al.ipOOS 



periodic boundary conditions. 

^3 



orm Cartesian g: 



rid 



with 

Our fiducial resolution 
features 256'^ zones, but some high resolution calcula- 
tions are done at 512^. We initialize the grid with uni- 
form fields. A background magnetic field is set along the 
a;-axis, with a magnitude Bq determined by the chosen 
value of /3 = pc3/[i?Q/87r], where p is the fluid density 
and Cs is its sound speed. The fluid velocity is initially 
zero. Kinetic energy is then steadily injected into the 
box, and the turbulence progressively grows. 

We drive the turbulence with a method similar to that 
of Lemaster fc Stone (2009). At each timestep, we gener- 
ate a velocity perturbation with a random amplitude in 
Fourier space, which is then added to the fluid's velocity. 
The Fourier components of this velocity perturbation are 
non-zero only for < fc < 4 x ^ (where L is the size of 
the simulation box), so that we are exclusively driving 
on large scales. We project each Fourier component per- 
pendicular to fc, so as to make it divergence-free. This 
avoids explicitly driving compressible modes. Finally, we 
ensure that the energy injection rate {E) is constant at 
each timestep; it is given roughly hy E ~ L^pSv^, where 
Sv is the magnitude of the turbulent velocity once the 
turbulence has saturated. 

After each random perturbation, the flelds evolve 
according to the conservative equations of isothermal 
MHD, either ideal or resistive, depending on the value 
of r]-. 



dp 
dt 



V • (pv) = 



We do not consider higher n resonances here. 



dpv 

Ik 



pvv 



BB 

in 



P 



B' 



= 



(6) 
(7) 
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dB 

~dt 



= V X X B) 



TjC 

in 



(8) 



In our simulations, these equations are implemented 
by evaluatin g the conservative fluxes w ith the HLLD Rie- 
mann solver ( Miyoshi fc Kusano|2005 ) and then evolving 
the fi elds using a Van Leer integrator ( Stone & Gardiner 
20091. The constrained- Transport (C'i'j algorithm en- 
sures that the magnetic field remains divergence-free. In 
the case of resistive MHD (77 ^ 0), the evolution of the 
fields is operator-split: at each timestep, we first evolve 
the fields according to the ideal equations, and then ap- 
ply the magnetic diffusion operator in a way that is con- 
sistent with CT. 

Our simulations are run with L = p = Cg = 1.0, but are 
actually scale-free. Our results can thus be transformed 
to any other physical parameters by scaling them with 
an appropriate combination of the new L, p and Cs- The 
behavior of the turbulence is controlled by two dimen- 
sionless numbers: the energy injection rate E (in units 
of pL'^C^) and the dimensionless magnetic field strength 
(e.g., (3 or va/cs). ^ 

Our choices of E and f3 produce both low sonic and 
low Alfvenic Mach numbers {Sv/cg < 1 , Sv/va < 1). 
Therefore, our simulations lie in the regime of weakly 
compressible Alfvenic turbulence. As the simulation pro- 
ceeds, the energy cascades from low wavenumber, where 
it is injected, to higher perpendicular wavenumber, thus 
forming an inertial range. The cascading energy is even- 
tually damped by physical and numerical dissipation at 
very high wavenumber. After some time (a few L/Sv), 
the energy dissipation rate equals the injection rate, and 
the amplitude of the turbulence saturates. Diagnostics 
of the power spectra in the saturated state agree rea- 
sonably well with the Goldreich-Sridhar predictions and 
with previous numerical simulations; e.g., our results are 



consistent with a k 



-10/3 



spectrum of the energy density 



and most of the energy is contained at perpendicular 
wavenumbers k± 3> k\\ , consistent with the critical bal- 
ance conjecture. 

3.2. The Particle Integrator 

Once the saturated state of the turbulence is reached 
(at t — AL/cs), we turn on the particle integrator. From 
then on, the particles and macroscopic fields are evolved 
simultaneously. The periodic boundaries of the simula- 
tion box also apply to the particles. We found that the 
driving of the turbulence could introduce spurious ac- 
celeration of the particles (see fjA] of the Appendix) . In 
order to avoid this we stop drivmg the turbulence after 
the particles have been injected, and let the turbulence 
progressively decay. The particle integrator usually runs 
for a physical time of l.OL/cs. Figurell] shows the evolu- 
tion of the turbulence during a typicalrun. The energy 
first builds up and saturates, and then decays as we inject 
the particles and turn off the driving. 

Our coUisionless charged particles evolve according to: 



clu 

Itt 



q 



mc 



u X B 



(9) 



where u is the particle's velocity and the remaining sym- 
bols have their standard meanings. In the particle inte- 
grator, the fields E and B are interpolated in space and 



0.12 




2 3 
Time ( x i/c,) 

Fig. 1. — Energy in tlie turbulent fluctuations versus time. Tlie 
solid and dashed lines correspond to the total kinetic energy and 
magnetic energy, respectively. The vertical dotted line marks the 
time at which the particle integration begins. After that time the 
turbulence decays, to avoid spurious particle heating due to the 
temporal driving (see Appendix rAll . 



time from their known values on the grid to the parti- 
cle's position. We obtain E from v and B on the grid 
by using: 



E 



1 



B 



47r 



V X B 



(10) 



The particle's charge-to-mass ratio (q/m in equation 
eq. [9]) is a free parameter. Our simulations simulta- 
neously evolve different types of particles, with charge- 
to-mass ratios spanning several orders of magnitude. 
To make the results easier to interpret, we convert 
the charge-to-mass ratio to a normalized gyrofrequency: 
r^o — qBo/mc. The actual gyrofrequency of a parti- 
cle is close to r^o, but may differ due to local turbulent 
fiuctuations in the B field. 

In our runs, particles all have positive charge. Since 
anisotropic Alfvenic turbulence consists of slow waves 
and linearly polarized Alfven waves, the sense of the gy- 
ration (due to the sign of the charge) should make no 
difference to the particle heating and acceleration. To 
confirm this, we compared simulations with negatively 
charged particles to those with positively charged parti- 
cles, and observed no significant differences. 

3.2.1. Integration method 

We first implemented the 4*''-order Runge-Kutta 
method and found that it was not suitable for this specific 
problem. We performed tests using a constant, uniform 
magnetic field (instead of the turbulent fields). Those 
tests showed a small but steady decrease in the gyrora- 
dius and the energy of the particles since the sign of the 
error is always in the same direction. With, e.g., 10 time 
steps per gyration, the energy decreased by 1 percent per 
gyration. This numerical bias is particularly dangerous 
given our goal of studying particle heating. 

We thus replaced the Runge-Kutta method by the fol- 
lowing implicit leap-frog numerical scheme, where posi- 
tions are defined at integer time and velocities at half 
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integer time: 
At 



(m"+1/2+u"-1/2) 



2+1/2 



At 



xB" (11) 
(12) 



where a superscript n indicates that the variable is de- 
fined at time level n. 



The specific implementation of equation (111 that we 
use is the Boris particle pusher ( Boris|[l970 1 MJnlike the 
Runge-Kutta method, this scheme is symmetric in time 
and symplectic, thus conserving energy and adiabatic in- 
variants to machine precision. Another advantage of this 
method is that the fields only need to be interpolated 
once per timestep, which makes it considerably faster 
than the 4*''-order Runge-Kutta method. 

We choose the particle timestep At so as to properly 
resolve the gyration. In order to also ensure that the 
particle timestep is smaller than the MHD timestep, we 
use: 

which produces at least 40 steps per gyration. Tests us- 
ing simplified field configurations confirmed the accuracy 
of this method. Some of our tests are described in i|B]of 
the Appendix. 

3.2.2. Interpolation 

The B and E fields are interpolated from their val- 
ues on the grid and at the discrete MHD timesteps to 
the particle's position. We use the Triangular Shaped 
Cloud (TS C) interpolation method ( [Hockney fc East^ 
wood|198T l in space and time. More details on the algo- 
rithm can be found in 5|C]of the Appendix. 

In ideal MHD, E'n = 0. Although this is indeed true on 
the grid, the interpolated E may have a non-zero com- 
ponent along the interpolated B. Because the particles' 
motion parallel to the magnetic field is unimpeded, this 
component, although small, can produce significant ac- 
celeration. In order to avoid non-physical parallel ener- 
gization, we compute E ■ B on the grid (this is zero in 
ideal MHD, but non-zero in resistive MHD) and inter- 
polate it to the particle's position. We then modify the 
component of E parallel to B: 



E 



{E-B)-E-B 



B 



(14) 



where the overline symbolizes interpolation and E is the 
new electric field. This ensures that: 



E-B = {E-B) (15) 
i.e., that the parallel component of E is correct. 

3.3. Measurement and Initialization of the Particles 

We initialize the particles uniformly over the entire 
computational domain. We then give them a parallel 
velocity W|| and a gyration velocity u± drawn from the 
distributions discussed below. 



Since the motion of a particle is in fact the superposi- 
tion of a rapid gyration and a slowly varying drift, cal- 
culating the gyration velocity u± requires knowledge of 
the drift velocity: u± = ||Mj^,tot — Ud\\. {u±.tot and 
are the total perpendicular velocity of the particle and 
its drift velocity, respectively.) We take into account the 
dominant E x B drift but neglect the other drifts (VB, 
curvature drift, polarization drift, ...), which are small 
for the particles considered here. For simplicity, we also 
neglect the E x B drift due to the resistive part of E,^ 
and adopt the following definition: 

u± = \\u±^tot - v±\\ (16) 

where Vj_ is the perpendicular velocity of the fiuid. 

We use two different kinds of initial velocity distribu- 
tions. For some of our calculations, our aim is to isolate 
the statistical effect of the turbulence on a certain type 
of particle. In this case, all of the particles have the same 
initial parallel velocity and magnetic moment: 

fa{x, u) ccS - U||^o) S (uj_ - iiqB{x)^ (17) 

where W|| q and /ig specify the initial particle properties. 
We give the particles the same magnetic moment, not the 
same gyration velocity, because analytic theory predicts 
that ^ should be conserved for particles with high fio 
{^^- With this initialization we can easily check whether 
the interaction with the turbulence indeed conserves /i. 
This test has been extensively performed while develop- 
ing our methods. 

By contrast, in other calculations, we study the inte- 
grated heating over a "realistic" distribution of particles. 
Here, we choose to assign the velocities according to an 
isotropic Maxwell-Boltzmann distribution: 



fMB{x,u) cx exp 



2c1 



(18) 



In all cases the particles are split into N groups, where 
in each group, particles all have the same physical prop- 
erties such as charge-to-mass ratio (i.e., fio) or initial 
velocity. Each of the N groups contains particles with 
different values of the physical property being studied in 
the particular simulation (e.g., fio or velocity). 

4. RESULTS 

In the following sections we describe the interaction 
between the turbulence and parti cles with different gy- 
rofrequencies fio and velocities (^4.11. We also com- 



pare simulations with and without explicit resistivity 
(i ]4.2| , simulati ons w ith different background magnetic 
field strengt hs (j :4.3), and simulations with different grid 
resolution (S4.4I. linally, we study the net effect of th e 
turbulence on a thermal distribution of particles ( ^4.5[ ). 
To help guide the reader's intuition, we note that our 
fiducial calculations (see Table jl]) have (3^1 and tur- 
bulence with outer-scale velocities 5v ^ O.Scg; the outer- 
scale eddy turnover time is ~ 2L/cs. We consider parti- 
cles with a range of velocities ^ 0. 1 — 10 . High velocity 
particles (3> Cs) transit the box many times in the course 
of the simulation; doing so, they may unphysically sam- 
ple a similar realization of the turbulent fluctuations. We 

^ This drift is proportional to »yc/47r, and is always small. 
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TABLE 1 

Properties of the Fiducial Calculation 



Parameter 


Value 


13 


1.0 


E 


O.lpL^cS a 


V 


(ideal MHD) 


Resolution 


256^ zones 


Integration time'' 


1.0 L/cs 


Number of particles'^ 


512 X 10^ 



^ This produces a sonic Mach number of ft: 0.46. Note 
that we stop driving the turb ulen ce once the particle inte- 
gration begins (see Fig. [l]& j ]3.2| . 

^ The time interval between the initialization of the par- 
ticles and the end of the simulation. 

The particles are initialized with a delta function in uy g 
and /ifl (cq. [17| ); we consider a range of ~ 10 — 
10* Cs/L. 



thus focus on lower velocity particles where the periodic 
box boundary conditions we employ are more reliable. 
This corresponds primarily to the velocities of protons 
or minor ions, rather than electrons, which have typical 
velocities ~ 40 Cs for ^ Tp. 

The range of gyrofrequencies we consider is ^ 10 — 
W^Cs/L. For a resolution of 256'^, the highest lin- 
ear frequency on the computational domain is ujmax — 
k\\^maxVA ~ 350 Cs/L where we have used the critical 
balance conjecture (Sj2| to estimate that for anisotropic 

Alfvenic turbulence, ^.max — kmaxkKf„ {kmax — Tr/Aa; 
is the highest wavenumber given a resolution Ax and 
krain — 27r/(L/4) is the wavenumber at which the tur- 
bulence is driven). Thus the particle gyrofrequencies we 
consider range from particles with gyrofrequencies com- 
parable to those of resolved turbulent fluctuations to par- 
ticles with 3> tOmax- From linear theory, the latter 
are expected to have primarily parallel heating because 
of the conservation of (S|2|. 



4.1. Fiducial Results 

The first calculations we describe are summarized in 
Tablejl] We consider an initial delta function of particles 
interacting with Alfvenic turbulence having Sv ~ 0.5cs in 
the presence of a mean magnetic field with P — iJ This 
interaction leads to both a stochastic diffusion in velocity 
space and a change in the mean energy of the particles; 
the magnitude of these changes depends on the cyclotron 
frequency f2o of the particles and their initial velocity. 
Figure |2] shows the standard deviation of u± and it|[, for 
U||,o = Cs and /io = c^/Bq (m_l,o — Cs), as a function of 
r^o, at the end of the integration (after l.OL/cs). Figure|3] 

Weakly compressible turbulence can nonlinear lv generate com- 
press ible fluctuations (fast waves in our case; e.g., |Cho &: Lazarian| 
|2003| ; this excitation is significantly weaker for smaller dv/Ca. Lo 
assess whether this could be energetically important in our case, 
we carried out simulations with a range of 5v/cs ~ 0.1 — 0.5 and 
found no significant differences relative to the results presented 
here, other than the expected scaling of the total particle heating 
and diffusion rates oc Sv'^. This rules out excitation of fast waves 
as an important source of heating in our simulations. 




10^ 10^ 

Qq ( X c JL) 

Fig. 2. — Standard deviation of the parallel My (solid line) and 
perpendicular u± (dashed line) velocities at the end of the simula- 
tion, as a function of the particles' gyrofrequency Qg; the particles 
initially have delta functions in velocity with M|| q = l.Ocs and 

/io = 1.0c^/-Bo {'>J-i_ — 1.0 Cs). The strong diffusion for low f!o 
particles is consistent with linear theory predictions for cyclotron 
resonance; a-(n||) 2> cr{u±^) for high f!o particles is consistent with 
linear theory predictions for /iV||i3 acceleration at the Landau res- 
onance (see i ]4.1| for details). The parameters of this calculation 
are summarized in Table 




10' 10^ 

( X c JL) 

Fig. 3. — Parallel (solid) and perpendicular (dashed) heating 
rates as a function of the particles' gyrofrequency Qq; the particles 
initially have delta functions in velocity with U|| q = 1.0 Cs and = 
1.0Cg/i?o («ix — 1.0 Cs). The heating rates are normalized by 
Eturb) the energy dissipation rate in the turbulence. The heating 
rates are averaged over the full l.OL/cs of the calculation. The 
parameters of this calculation are summarized in Table ^ 



shows the parallel and perpendicular heating rates for the 
same calculations (l^||,± = s Epartzde. s'^^f ■ To 
make the heating rates easier to interpret, we normalize 
them by the average energy dissipation rate in the de- 
caying turbulence. Because our calculations are for test 
particles, the heating rate can be scaled to any ptest/p, 
i.e., to any test particle density. 

The dispersion (Fig. [2| and heating (Fig. |3| pro- 
duced by the turbulence depend on the gyrofrequency 
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of the particles. Particles with low gyrofrequency fip ^ 
300 Cs/i are the most strongly heated, with most of the 
heating being in the perpendicular direction (i.e., u± in- 
creases rather than uy). The heating also depends on 
the gyrofrequency of the particle. By contrast, particles 
with high gyrofrequency f2o ^ SOOc^/i have a signifi- 
cantly lower heating rate that is independent of ^Iq; in 
addition, most of the dispersion/heating is in the direc- 
tion of the magnetic field. We distinguish between the 
two classes of particles in Figures |2]& [3] by characterizing 
the low r^o particles as "cyclotron-resonant" and the high 
fig particles as "Landau-resonant" (the reasons for these 
particular identifications will become clear shortly). 
The significant heating of low fio particles in Figures 

f&[3]is consistent with the analytic expectations for cy- 
otron resonance (^. To confirm this, we note that the 
range of Alfven frequencies present in the simulation is 
~ 4 — 350 Cs/L. This is comparable to the range of 
in Figures |2] & [3] over which there is strong heating and 
diffusion. The heating is also primarily perpendicular 
heating, as expected analytically. The negative parallel 
heating for the cyclotron-resonant particles in Figure [3] 
is a consequence of pitch-angle scattering, as we demon- 
strate in more detail below. 

Figures |2] & [3] show that particles with high gyrofre- 
quency (rip ^ 10^ Cs/L) are significantly less affected by 
the turbulence than the low fip cyclotron-resonant par- 
ticles. This is particularly true for the perpendicular 
velocity u±: cr(u±) <C o'(u||) at high Hq. We interpret 
this low dispersion as a nonlinear analogue of the absence 
of perpendicular heating in linear theory (|2|, which is 
ultimately due to the adiabatic invariance of the mag- 
netic moment. In fact, a{u±) at high fio in Figure [2] is 
only modestly larger than the initial dispersion in u±; 
because we initialize the particles with a fixed fj,, there is 
an initial dispersion in u± due to differences in B in the 
turbulent plasma. 

To assess the changes in n explicitly. Figure 4] shows 
the standard deviation of ^ and the change in the mean 
of fj, for the same calculations as in Figures |2] & |3] At 
the end of the calculation, after 1.0 L/cg, the mean ^ of 
the Landau- resonant (large fig) particles has changed by 
< 1%; the dispersion introduced by the interaction with 
the turbulence is somewhat larger ~ 10%. These small 
changes are qualitatively consistent with the prediction 
of linear theory that high gyrofrequency particles evolve 
adiabatically and thus conserve their magnetic moment. 
Nonetheless, the small changes in /i in Figure |4] appear 
to be real and represent a quantitative difference rela- 
tive to linear theory. We have carried out a number of 
tests to confirm that the changes in /i shown in Figure 
[4] are not numerical; e.g., the results are independent of 
the timestep in the particle integrator and of the Courant 
number in the integration of the turbulence. In addition, 
tT(/i) cx 5v, consistent with how the parallel diffusion de- 
pends on the amplitude of the turbulence; this argues 
against violation of conservation due to finite ampli- 
tude waves, as has been proposed in other contexts (e.g., 
Johnson & Cheng||2001 1. One possible explanation for 




10^ 10^ 

( X c JL) 

Fig. 4. — Standard deviation of the magnetic moment at the 
end of the simulation (dashed Une) and the change in its mean 
(sohd line), as a function of the particles' gyrofrequency f2o; the 
particles initially have delta functions in velocity with un q = 1.0 

and fiQ = 1.0c^/_Bo ("±,0 — l-Ocs). The large changes in fi at low 
Co are due to cyclotron resonance. Linear theory predicts that 
is conserved for high Hp , which is not fully consistent with the 
numerical results; see §4.1| for an interpretation. The parameters 
of this simulation are summarized in Table [l] 




the changes in u± is that the delta function resonances in 
linear theory (eq. [l]) are significantly broadened in the 
nonlinear turbulence (as we show in more detail shortly) . 
Recall from equation Q that the change in u±_ due to 
interaction with a single Fourier component of the tur- 



FlG. 5. — Heating rates in the parallel (solid line) and perpen- 
dicular (dashed line) directions as a function of the initial parallel 
velocity M||_0) foi" Landau-resonant particles {Qq = 4.0 X 10^ Cs/L). 
Although varies, all particles initially have the same magnetic 
moment, fio = I.Oc^/Sq (and thus uj_ q ~ l.Ocs). Linear theory 
predicts that only particles with M|| g = ~ O.Scs (the Landau 

resonance) should exchange energy with the waves. By contrast, 
particles with all M||_o, and particularly those with q < Cs, inter- 
act strongly with the turbulence. The heating rates are normalized 
by the total initial energy of the particle Etot,i = + ^^'^ 
are averaged over the whole integration. The properties of this 
simulation are summarized in Tabled 

bulence is du±/dt oc [lu — Iii linear theory this 

vanishes because non-zero energy exchange with a fluc- 
tuation requires lo = fc||U||; this restriction is relaxed in 
strong Alfvenic turbulence (e.g., Chandran|2000 | so that 
small changes in uj_ may be possible. 

As described in ^ analytic theory predicts that high 
r^o particles interact with anisotropic Alfvenic turbulence 
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Landau-resonant Cyclotron-resonant 




Fig. 6. — Scatter plots in velocity space (My, uj_) for Landau-resonant particles (Cq = 2.0 X 10*Cs/L; left panels) and cyclotron-resonant 
particles {Qq = 60 c^/L; right panels) at the end of the simulation. The initial distribution function is a delta function. Landau-resonant 
particles primarily diffuse in ny while cyclotron resonant particles undergo pitch angle scattering which tends to isotropize the distribution 
function. In the upper panels, the particles initially have low velocities, u^fi = U|| q = 0.3cs. In the lower panels, the particles initially 
have high velocities, uj_ q = q = 3.0cs. The remaining parameters are summarized in TablelTl 



primarily via /iV||i3 acceleration by the slow magne- 
tosonic modes. Anisotropic slow waves have t he disper- 



sion relation w ~ Cs?;A|fc|||(Cs -I- vj^) 



-1/2 



(e.g., 
- 0.8 



Lithwick 

ku Ic, for 



& Goldreich 2001), which reduces to w 
/j = i. ihus tor our fiducial calculation, the condition 
for high flo Landau resonance (n = in eq. [l]) becomes: 

M|| ~ ±0.8c,. (19) 

Note that the resonant condition does not depend on k, 
but only on uy: the resonant particles simultaneously 
interact with all modes, but only a single parallel ve- 
locity is picked out in linear theory. Figure [5] shows 
the parallel and perpendicular heating as a function of 
the initial parallel velocity U||^o for our fiducial calcula- 
tion, taking a fixed /iq = 1.0c^/Bq and a sufficiently 
high r^o = 4.0 X 10^ Cg/L that there is no cyclotron res- 
onance. Contrary to linear theory, particles with a wide 
range of velocities 0.1 — 10 receive significant energy; 
in particular, for U|| g ^ O.Scs (the resonant velocity in 
linear theory), the heating is relatively independent of 
M|| g.® The wide range of velocities at which the particles 
couple to the turbulence is qualitatively consistent with 
the idea that the linear resonances (eq. [l]) are highly 
broadened in Alfvenic turbulence because the nonlinear 
decorrelation time is comparable to or shorter th an the 
linear peri od of the Alfvcn and slow waves (e.g., Chan- 
|dran"2000). We defer a more detailed test of resonance 
Droadening models to future work. Although the linear 

* Note that the results for Figures jSj and js] were for g = Cs, 
for which the heating is particularly small. 



resonance no longer manifests itself as a delta function. 
Figure [5] shows that at nearly all velocities, the heat- 
ing of high rig particles is primarily parallel to the local 
magnetic field, consistent with linear theory predictions. 

Figure [6] shows the positions of the individual parti- 
cles in velocity space at the end of the fiducial simula- 
tion, for both high fig Landau-resonant particles (left 
column) and low f2g cyclotron-resonant particles (right 
column). The top row is for low velocity particles with 
U|| g ~ uj_fl ~ 0.3 while the bottom row is for high 
velocity particles with U|| g ~ w_L,g ~ Sc^. For Landau- 
resonant particles, the dispersion occurs almost exclu- 
sively in the parallel direction, as we have previously 
seen. For the cyclotron-resonant particles, the diffu- 
sion in velocity space depends somewhat on the velocity 
of the particles. For uy g < va (top right), the par- 
ticles are preferentially dispersed in the perpendicular 
direction due to the cyclotron resonance. However for 
U||^g > VA (bottom right), the velocity diffusion is domi- 
nated by pitch angle scattering: the velocity distribution 
is isotropized more rapidly than the total energy changes. 
This difference in the results of cyclotron resonance for 
sub and super-Alfvenic p articles is a well-known result 
of quasihnear theory (e.g.. Miller fc Roberts||l995 1. 



4.2. Resistive Simulations 

In the ideal MHD simulations described in the previous 
sub-section, i^y = because of our constrained transport 
algorithm; this is also preserved to machine accuracy by 
our interpolation methods (§31). Numerical reconnection 
is nonetheless present on small scales. To understand the 
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Fig. 7. — Standard deviation of the parallel My (top panel) and 
perpendicular nx (bottom panel) velocities at the end of the simu- 
lation, as a function of the particles' gyrofrequency Oq- The curves 
are plotted for different resistivities which are converted to mag- 
netic Reynolds numbers, as computed on the outer scale of the tur- 
bulence. The solid line is for r] = (ideal MHD), while the dashed, 
dot-dashed and dotted lines have increasing resistivity. The pres- 
ence of 7^ produces an that leads to significant parallel 
heating and diffusion at high flo , but the properties of the low f!o 
cyclotron-resonant particles do not depend significantly on 77. The 
particles initially have delta functions in velocity with M|| q = 1.0 Cs 

and /XQ = l.Oc^/Bg (lix — l.Ocs). Aside from rj, the parameters 
of this simulation are the same as in Table [l] 

effects of reconnection and small-scale current sheets on 
our results in a more controlled manner, we have carried 
out a number of simulations with finite resistivity rj.^ It 
is important to stress that the physics of reconnection is 
not adequately represented by a spatially and temporally 
constant value of 77. Thus our calculations with finite 77 
should not be interpreted as physical. Rather, they allow 
us to isolate the potential importance of reconnection for 
some of our results and assess whether the interpretations 
given in the previous sub-section are correct. 

We parameterize our chosen values of rj using the mag- 
netic Reynolds number of the turbulence on the outer 

^ We keep the energy injection rate E constant, which leads 
to a slightly lower Mach number for larger rj, due to the higher 
dissipation. 



scale: Rm = L5v /[rjc/^ir]. We chose values of 77 so that 
Rm ~ 1 on the grid scale, i.e., so that small-scale current 
sheets are marginally resolved; this implies 



A a; Svh 



1 



(20) 



where Ax is the size of a cell and k^ax — tt/Ax. We 
estimate Sv^^^^ assuming a Kolmogorov cascade Sv^ ~ 
5v{kL)~^/^ . For our fiducial resolution of 256^, Rm ~ 1 
on the grid scale thus corresponds to Rm ^ 1500 at the 
outer scale of the turbulence. 

Figure [7] compares the velocity dispersion in the paral- 
lel and perpendicular directions at the end of simulations 
with different values of rj, for particles with My — 1.0 Cs 
and /iQ = 1.0 c^/i?. The most striking result is the very 
strong parallel dispersion for particles with high gyrofre- 
quency. This is a consequence of the non-zero parallel 
electric field in resistive MHD: E\\ = ^(V x As 
a result, particles are freely accelerated in the parallel 
direction by the qE\\ force (see eq. [sj). The paral- 
lel acceleration thus depends strongly on the particle's 
charge-to-mass ratio q/m - or equivalentlvon its f2o- 
This dependence is clearly visible in Figure]?] 

The perpendicular velocity u±^ is much less affected 
by resistivity, since the additional perpendicular electric 
field Eji_ — j^(V X B)x^ averages to zero over a gyration. 
The strongest effect evident in Figure |7] is that, for fio ^ 
lOOcs/i, the perpendicular diffusion becomes weaker as 
the resistivity increases. The most probable explanation 
is that resistive dissipation preferentially damps the high 
k modes, which are precisely the modes that resonate 
with particles having f2o ~ lOOc^/i. 

Overall, our simulations with an explicit resistivity 
demonstrate that, in resistive MHD, current sheets lead 
to parallel heating by the non-zero E\\ . This is in addition 
to the parallel heating by iiV\\B forces and the perpen- 
dicular heating and pitch-angle scattering by cyclotron- 
resonant waves highlighted in the previous sub-section. 

4.3. Variations in (3 

We carried out a number of calculations at different /3 
to asses s whether the particle heating physics identified 
in §4.1| depends significantly on (3. Because we are fo- 
cusing on weakly compressible Alfvenic turbulence, we 
kept the Alfvenic Mach number Sv/va at the saturation 
of the turbulence constant when varying /3 (by appropri- 
ately varying E). There were no significant changes in 
any of our results as a function of f3. To illustrate one 
example. Figure [8] shows how the standard deviation of 
u± at the end of the simulation depends on the back- 
ground field. In this calculation, the particles initially 
have jiQ = 1.0c^/-Bo and uy q = 0, the latter to avoid 
any Doppler shifts in the cyclotron resonance. 

For cyclotron resonant particles, the decreasing mag- 
nitude of a{u_L) for higher /3 in Figure[8]is due to the fact 
that, to keep Sv/va constant, the energy in the turbulent 
fiuctuations is lower for high /?. The cyclotron resonance 
clearly shifts towards higher finas (3 decreases. This can 
be understand from equation ([5| and the assumption of 

We found similar results in ideal MHD, when we did not ex- 
plicitly constrain to be zero when interpolating from the MHD 
grid to the particle's position (see ^3.2.2[ l. 
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Fig. 8. — Standard deviation of the perpendicular velocity uj_ 
at the end of the simulation as a function of the particles' gy- 
rofrequency VIq, for ny q = (to minimize Doppler shift) and 
u± Q ~ l.Ocs. The different curves correspond to different val- 
ues of 13. Linear theory predicts that cyclotron resonance with 
Alfven waves requires ~ Qq, consistent with the increasing 

heating/dispersion at high Oq for decreasing (3 (i.e., increasing v^). 
E is adjusted so that all calculations have approximately the same 
Alfvenic Mach number. 

Alfvenic turbulence: resonance requires ~ ilo and 

thus the resonant particles should have fig oc This 
scaling is reasonably consistent with the shift in Figure 



4.4. The Effects of Numerical Resolution 



Finite numerical resolution limits the conclusions we 
can draw for several reasons, even given our restriction 
to scales where MHD is valid. In particular, higher reso- 
lution changes the properties of the turbulence, since we 
can resolve higher k modes. For the same reason higher 
resolution increases the range of linear and nonlinear fre- 
quencies found in the turbulence, and thus the range of 
particles that can be cyclotron resonant. 

Figure [9] shows the parallel and perpendicular heating 
rates for different resolutions. Going to higher resolution 
does not significantly affect the parallel heating rate at 
high flf). Th i s is re asonably consistent with linear the- 
ory: 



Barnes 



([1966|) showed that the linear coUisionless 
damping of the slow magnetosonic mode is independent 
of the magnitude of k (although it depends on its direc- 
tion). Thus increasing the range of k contained in the 
domain should not significantly change the heating rate, 
since all scales contribute significantly. By contrast, Fig- 
ure[9]shows that the range of gyrofrequencies fig that can 
be cyclotron resonant increases with resolution: higher 
frequency particles are cyclotron-resonant at higher reso- 
lution because the range of resolved k, and thus the range 
of resolved frequencies, increases at higher resolution. 

Quantifying the range of f2o that are cyclotron- 
resonant at a given resolution is critical for understand- 
ing the implications of our results for solar and astro- 
physical problems (35.21. Towards this end we define a 



maximum cyclotron frequency flmax'- this is the value 
of r^o at which Afi/fi, a{uj_), or cr(/i) is larger than its 
flo oo value by a factor of 2. The exact definition of 
^max is somewhat arbitrary, but our definition captures 
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Fig. 9. — Parallel (upper panel) and perpendicular (lower panel) 
heating rates, as a function of the particles' gyrofrequency f!o> 
for M|| Q = l.Ocs and u^ q ~ l.Ocs. The different curves corre- 
spond to different resolutions. With increasing resolution, higher 
frequencies are present on the computational domain, which in- 
creases the range in f!o over which there is perpendicular heating 
due to the cyclotron resonance. By contrast, the parallel heating at 
high Co is relatively independent of resolution. The heating rates 
are averaged over the whole integration (i.e., for 1.0 L/cs), and 
are normalized by Bturbi the energy dissipation rate in the turbu- 
lence. Apart from resolution, the properties of these calculations 
are summarized in Table (T) 



the key result seen in many of the plots in this paper 
(e.g.. Fig. [2]&l4|: the dispersion/heating in /i and u±^ 
are small athign fig (for Landau-resonant particles) and 
then increase rapidly for smaller fig as f2o becomes com- 
parable to the frequency of Alfven waves on the compu- 
tational domain. By considering three different physical 
quantities (A/i//i, (t(u_l), and cr(/i)) in defining 0.max we 
hope to bracket some of the uncertainty in our estimate 

of ^rnax- 

Figure [lO] shows these three different definitions of 
^max as a function of h 
calculation with /3 = 1, wy q 
There is a clear tend of inCrGclsing ^max with ^maa:- 

Also 

shown in Figure [TO] are two theoretical predictions for 
the maximum Alfven wave frequency as a function of 
resolution. The first (dotted line) assumes an isotropic 



max = tt/Ax, for our fiducial 
Mil n — 1.0 Cs, and u_l,o — 1.0 Cg. 
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Fig. 10. — Maximum gyrofrequency Slmai for cyclotron-resonant 
particles as a function of resolution (colored symbols); the resolu- 
tion is quantified by the maximum wavenumber resolved in the 
simulation kmax = tt/Ax. Also shown is the maximum frequency 
of Alfven waves as a function of resolution in an isotropic cascade 
(dotted line) and for anisotropic Alfvenic turbulence (solid line; 

~ k^axk^^f^). The numerical results favor the anisotropic 
turbulence models. Qmax is the value of the particle gyrofrequency 
Qq at which Afi/fi, aOuj), or cr(fi) is larger than its Qg — > oo value 
by a factor ^ 2 (see j]4.4|for details). 
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Fig. 11. — Distribution functions vs. total velocity u = 

^M'^ -I- My , for cyclotron-resonant particles with = 250 Cs/L. 

The distribution function is initially thermal with the same sound 
speed as the fluid {solid line); the test particles thus represent pro- 
tons. The dashed line is the distribution function at the end of 
the integration (after l.OL/cs), while the dotted line is a thermal 
distribution with the same average energy as the final distribution. 
The cyclotron resonance leads to modest acceleration of high en- 
ergy particles. The properties of the simulation are the same as in 
Tablenl except for the resolution, which is 512'^. 



cascade in which case the maximum Alfven frequency is 
— k^ax^A- The second estimate (sohd hne) takes into 
account the anisotropy of Alfvenic turbulence; critical 
balance implies that the maximum Alfven frequency is 



k 



,2/3 , 



1/3 



va, where in Figure 



10 



we take 



kmin — 27r/(L/4) sincc we drive modes with wavelengths 
between L/A and L (see S3.ll. 

The estimates oi^max in 1 igure 10 are reasonably con- 
sistent with ^max — '^k\\,maxVA, which is comparable 
to the maximum Alfven wave frequency in anisotropic 
Alfvenic turbulence. Because the discrete resonances of 
linear theory are "broadened" in strong MHD turbu- 
lence (Fig. p|, the maximum gyrofrequency of particles 
that feel the cyclotron resonance should be somewhat 
larger than maximum Alfven wave frequency. This is 
true for the anisotropic estimate of ^max in Figure |10| 
but not the isotropic estimate. The non-trivial point 
demonstrated by Figure 10 is that (as predicted by lin- 
ear theory) it is the smallest parallel length-scale in the 
turbulent fluctuations that determines the efficacy of cy- 
clotron resonance, and thus perpendicular heating. This 
fact will be very important when we apply our results to 
space physics and astrophysics problems in §5.2[ 



4.5. Heating of a Thermal Distribution 

In the previous sections we have focused on the diffu- 
sion and heating of particles that all start with approxi- 
mately the same velocity, in order to isolate the physics 
of test particles interacting with strong Alfvenic turbu- 
lence. In order to obtain a more realistic estimate of the 
total heating rate in a turbulent plasma, we now con- 
sider a n is otropic thermal distribution of particles {Jmb 
in eq. 18 ). These calculations have (3 = \ and the test 



particles have a thermal speed equal to that of the fluid; 
thus the thermal test particles represent protons or other 




10-' 10^ 
( X c JL) 

Fig. 12. — Parallel (solid line) and perpendicular (dashed line) 
heating rates, as a function of the particles' gyrofrequency VIq. The 
heating rates are averaged over the last half of the simulation (i.e., 
for 0.5 L/cs), and are normalized by Eturh^ the energy dissipa- 
tion rate in the turbulence. The distribution function is initially 
thermal with the same sound speed as the fluid; the test particles 
thus represent protons. The Landau-resonant particles with high 
Co receive primarily parallel heating while the cyclotron-resonant 
particles with lower Ho are largely heated in the direction perpen- 
dicular to the local magnetic field. 

particles with a similar thermal speed. 

Figure [TT] shows the evolution of the velocity distribu- 
tion as the simulation proceeds, for cyclotron-resonant 
particles (fig = 250cs/i). The distribution function de- 
velops a modest non-thermal tail, although most of the 
energy remains in the thermal population of particles. 
Figure [12] shows the particles' parallel and perpendicular 
heating rates, divided by the average energy dissipation 
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rate in the turbulence, as a function of the particles' gy- 
rofrequency fio- As before, these results can be rescaled 
by Ptest/p to apply them to any specific population of 
test particles. For high gyrofrequency (> SOOc^/L), 
Figure [12] shows that the heating is primarily parallel 
to the magnetic field; the total perpendicular heating 
rate fiuctuates around ~ —0.025 i?t,jrb — ~-E'||/20, suf- 
ficiently small that it is hard to read off of the Figure. 
As we have argued before, this parallel heating is due to 
the stochastic ^VB forces created by slow magnetosonic 
waves. Note that the magnitude of the parallel heating 
rate is ~ 0.5 Eturb- This is sufficiently strong that, for the 
parameters chosen here (/3 = 1), a significant fraction of 
the energy in the slow wave cascade would be dissipated 
on large scales, rather than cascading to small scales. 
For lower fig (10 — 300 Cg/i), Figure 12 shows that the 



net heating of the thermal distribution function is pri- 
marily perpendicular to the magnetic field. The heating 
is very strong, ~ SEturb (for ptest — p) and ~ 6 times 
the heating rate in the high fio limit. This strong per- 
pendicular heating is, we have argued, due to cyclotron 
resonance. For ptest — P, the energy gained by the parti- 
cles significantly exceeds the energy available in the tur- 
bulence. This clearly demonstrates that the damping is 
sufficiently strong that the test particle approximation 
breaks down, if in fact the turbulent fluctuations reach 
w ~ r^o- We return to this issue in f|5j 

5. SUMMARY AND IMPLICATIONS 

We have carried out a detailed study of the heating of 
test particles by weakly compressible MHD turbulence. 
Our calculations integrate the equations of motion for 
up to ~ 10^ particles in "real time" as the turbulence it- 
self evolves. The heating and acceleration of particles by 
MHD turbulence plays a central role in theoretical mod- 
els of many heliosphe ric and astrophysical ph enomena, 
including solar flares (Miller & Roberts||1995||, cosmic- 



ray sc attering and confinement in the Galaxy (iChandran 



2000|, the hea ting and acceleration of the solar wind 



( [Cranmer fc van Ballegooijen 2005), and the radiation 
produced by some accretion disks onto compact objects 
(Quataert & Gruzinov 19991. Because of the many po- 
tential applications of this work, we have focused on the 
basic physics of particles interacting with MHD turbu- 
lence, rather than on any of these specific applications. 
In this section we first summarize our primary results 
and then briefiy discuss their implications for astrophys- 
ical systems, in particular the solar wind. 

Our turbulence is driven subsonically, sub- Alfv enically 
and with a divergence free velocity field ([3.11. Such 
weakly compressible MHD turbulence consists of nonlin- 
early interacting Alfven and slow magnetosonic waves, 
with most of the energy cascading to sma ll scales per- 
pendicular to the local magnetic field (e.g., Goldreich & 
Sridhar| 1995|. There is an extensive body of work study- 
ing the heating of test particles using (quasi)linear the- 
ory, which predicts th at energy e xchange happens at dis- 
crete resonances (e.g., Stix||1992|an d references therein); 
we review these predictions in ^ Our calculations re- 
lax many of the simplifying assumptions made in linear 
theory. Specifically, rather than studying the interac- 
tion between particles and a superposition of long-lived 
linear waves, our particles interact with strong MHD tur- 



bulence, i.e., with turbulence in which the timescale for 
nonlinear interactions to transfer energy to smaller scales 
is comparable to, or shorter than, the linear period of the 
waves. We have also carried out both ideal and resistive- 
MHD simulations in order to isolate the importance of 
dissipation in current sheets, rather than wave-particle 
resonances, for particle heating. 

5.1. Summary 

One of the striking conclusions of our work is that 
many - although not all - of the results we find for test 
particle diffusion and heating in fully developed MHD 
turbulence are qualitatively consistent with the predic- 
tions of quasilinear theory (we defer a more quantitative 
comparison between our numerical results and quasilin- 
ear theory to future work). More specifically, our pri- 
mary results can be summarized as follows: 

• How particles are heated depends sensitively on 
whether the gyrofrequency of the particle is compa- 
rable to the frequency of a turbulent fluctuation lj that 
is resolved on the computational domain. 

• Particles with ^ lo undergo strong perpendicular 
heating and pitch angle scattering, qualitatively consis- 
tent with linear theory predictions for cyclotron reso- 
nance with the linearly polarized Alfven wave s present 



121. In \ b.2 



in MHD turbulence (e.g.. Fig. |2| 

we discuss the implications of tms perpendicular heat- 
ing for measurements of proton and ion temperature 
anisotropics in the solar corona and solar wind. 

• Particles with fip ^ undergo strong parallel heating, 
qualitatively consistent with linear theory predictions for 
the heating produced by the /iVyi? forces ("transit time 
damping") associated with the slow magnetosonic waves 
in MHD turbulence (e.g.. Fig. |2] [S] [6] & [l2]) . 

• In contrast to the predictions of linear theory, we 
find that discrete resonances do not dominate the energy 
transfer between particles and waves. Instead, energy 
transfer occurs for a wide range of particles, even those 
that would be quite "off-resonance" according to linear 
theory (Fig. |5|. This is broadly consistent with models 
in which the rapid decorrelation of waves in anisotropic 
MHD turbu lence leads to sig nificant "resonance broad- 
ening" (e.g., |Chandra^|2000[ ). 

• Linear theory predicts that particles with fio ^ un- 
dergo purely parallel heating because of the adiabatic 
invariance of the magnetic moment p, cx u\lB in the 
presence of low-frequency electromagnetic fluctuations. 
Although we do find that most of the heating is paral- 
lel to the local field for particles with fio ^ w (Fig. |3] 
& |12|) we also see small, but non-zero, changes in u±^ 



and p even for high cyclotron-frequency particles (Fig. 
[4]). The physical origin of this perpendicular heating is 
not fully clear, but we speculate th at it may be a con- 
sequence of resonance broadening ([4.11. We find that 
the perpendicular heating and diffusion rates are ~ 10% 
of the parallel heating and diffusion rates for fio 3> 
(e.g.. Fig. [5] fc [l2|. Although this is unlikely to change 
any of the general conclusions about perpendicular ver- 
sus parallel heating derived from linear theory, it is an 
interesting modification to the physical picture of how 
particles interact with anisotropic Alfvenic turbulence. 
• All of the results summarized above are present in both 
ideal MHD simulations and in simulations that include 
an explicit resistivity to dissipate small-scale fluctuations 
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in the turbulence. A finite resistivity also generates a 
non-zero that produces strong parallel heating and 
diffusion of particles (Fig. [?]). It is important to stress 
that the physics of reconncction is not adequately rep- 
resented by the spatially and temporally constant resis- 
tivity present in our calculations. Indeed, particle-in-cell 
and Hall MHD simulations of reconnection show that 
particle heating in current sheets is far more subtle than 
can be captured by resistive-MHD simulations (D rake] 
|et al.||2009b|a| ). Further work is thus required to un- 
derstand the particle heating and acceleration produced 
by current sheets that self-consistently develop on small 
scales in a turbulent plasma. 

• For the particular case in which we initialize a thermal 
distribution of test particles having the same sound speed 
as the fluid in our turbulence calculations, the test par- 
ticles approximate the interaction between protons and 
the turbulent fluctuations. For (3^1, Figure [T2| shows 
that the heating rate of the test particles is ~ •iEturb for 
cyclotron-resonant particles with low CIq and ~ O.bEturb 
for Landau-resonant particles with high fio (where Eturb 
is the energy dissipation rate in the turbulence). This 
highlights that in a low-coUisionality plasma, the wave- 
particle interactions are strong enough to damp a large 
fraction of the turbulent energy; the test particle approx- 
imation is thus a poor one. 

In all of our calculations, the test particles have ve- 
locities within a factor of '-^ 10 of the sound speed of 
the fluid Cs- This corresponds primarily to the velocities 
of protons or minor ions, rather than electrons, which 
have typical velocities ~ 40 Cg for ^ Tp. High velocity 
particles (3> c^) transit the computational domain many 
times in the course of the simulation; doing so, they may 
repeatedly sample very similar realizations of the tur- 
bulent fluctuations. We believe that tests with larger 
computational domains are required in order to reliably 
study the turbulent diffusion and heating of high veloc- 
ity particles. This limitation has two significant implica- 
tions. First, our calculations do not at this point predict 
the relative turbulent heating of electrons and protons, 
which would be of considerab le interest for both solar and 
astrophysical problem s (e.g., Quataert fc Gruzinov|1999 



Cranmer et al. [2009 1 . Second, our calculations are lim- 
ited in their ability to predict the acceleration of particles 
to energies well above the thermal energy of the plasma. 
It is nonetheless important to highlight that when we 
initialize a thermal distribution of test protons that are 
cyclotron resonant with the turbulent fluctuations, the 
distribution function naturally evolves a non-thermal tail 
(Fig 
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Previous work o n test particle heating in MHD turbu- 
lence simulations (Dmitruk et al. 20041 concluded that 
current sheets produce significant parallel heating of elec- 
trons and that ions are preferentially heated in the per- 
pendicular direction. Our results are broadly consis- 
tent with these conclusions give n Dmitruck et al's as- 
sumed particle gyrofrequencies. Dmitruk et al. (20041 



give an involved physical interpretation of the perpen 
dicular heating in their calculations; we suspect that it 
is simply due to cyclotron resonance, as we have found in 
our calculations.^""^ Perhaps most importantly, the per- 

Although [Dmitruk et ar| | |2004[ | studied test particle heating 



pendicular io n heating found both here and in Dmitruk 
et al. (20041 is not, we believe, applicable to the solar 
wind, contrary to the claims ma de by Dmitruck et al . 
In particular, as we now explain, Dmitruk et al. (20041 
did not consider the limitations imposed by finite numer- 
ical resolution when claiming that their results could be 
directly applied to the solar wind. 

5.2. Implications 

Broadly speaking, our results suggest that many of 
the conclusions drawn from quasilinear theory about 
the heating and acceleration of particles by anisotropic 
Alfvenic turbulence are likely to be qualitatively cor- 
rect. However, many of the quantitative results may 
change given, e.g., the lack of the discrete resonances 
that strongly shape the predictions of quasilinear the- 
ory (Fig. Assessing in detail the implications of our 
results for heating by MHD turbulence in solar and astro- 
physical environments will require additional work. Here 
we take the near-Earth solar wind as an example to illus- 
trate the implications of our results and the limitations 
due to finite numerical resolution. 

A number of observations indicate that the solar 
wind undergoes spatially extended heating. For exam- 
ple, in situ measurements from satellites such as He- 
lios & Ulysses show that electrons an d protons have 
non-a diabatic temperature profiles (e.g., Cranmer et al. 
20091. Heating by anisotropic Alfvenic turbulence is 



adiabiticity (e.g., 


Matthaeus et al.|1999 


Cranmer & van 


Ballegooijen |2003 


p. In situ measurements of the solar 



wind at i AU snow tiiat the proton distribution lunc- 
tion is on average anisotropic with respect to the local 
magnetic field: ~ 0.9 Ty ( ^Bale et al.,,2009) , although 
the sign of the anisotropy depends on the wind speed, 
with Tj. > Tn for v^ind > 600 km s'^ and Tj. < Tn for 



jjumd ^ 600 km s ^ ( [Kasper et al.||2002 Helhnger et al. 

; clear how to understand these mea- 



2006|. It is not yet 
sureo temperature anisotropies in the context of heating 
by anisotropic Alfvenic turbulence. 

Typical physical parameters for the slow solar wind^^ 
at - 1 AU near Earth are B ~ lO"'' G, no - 20 cm'^, 
Tp~Te~ 1.5 X 10^ K, VA ^ 50 km s-\ - 35 km s~\ 
z^wind — 460 km s""", and (i ~ 0.4 (e.g., Celnikier et al. 
19871); the proton Larmor radius and gyrofrequency are 



thus p ~ 3 X 10^ cm and flp ~ 0.15 Hz, respectively, 
while the electron Larmor radius and gyrofrequer icv are 
e — 10^ cm and fie — 300 Hz, respectively. Howes 



et al. ( 2008a I reviewed a variety of observational diagnos 



tics of the outer scale of the turbulence in the solar wind 



and estimated that fc„ 



10 



10-11 



cm. The cyclotron 



frequency in our calculations is expressed in units oicg/L 
where the size of our box L ^ in /kmin is also of order the 
outer-scale of the turbulence Expressed in these 

units, the proton and electron cyclotron frequencies in 

in a static snapshot of MHD turbulence (i.e., uj = 0) their parti- 
cles could nonetheless undergo cyclotron resonance via the parallel 
Dopplcr shift in eq. Hv. 

Our simulations nave roughly the same energy flux traveling 
in both directions along the local magnetic field. This is typically 
not true at ~ 1 AU in the solar wind, particularly in the fast wind 
fMarsch 2006). For this reason, we compare to measurements in 
the slow solar wind, where the assumption of balanced turbulence 
is more consistent with the measurements. 
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the solar wind are ~ 10'' ^ Cg/L and 



10 



7-8 . 



s/L, 



when the cascade reaches k i 



is fir; 



0.02 



spectively. The minimum value of the proton cyclotron 
frequency in the solar wind is thus comparable to the 
maximum value of the cyclotron frequency of particles 
that we have simulated (e.g., Fig 12 1. The fundamen- 
tal reason for this is that, in nearly all heliospheric and 
astrophysical plasmas, the ratio of the proton cyclotron 
frequency to the outer-scale frequency of MHD turbu- 
lence is much larger than the dynamic range that can be 
simulated with current computational resources. 

A naive application of our results to the near-Earth 
solar wind, using the parameters above, suggests that 
all particle heating by MHD turbulence corresponds to 
Landau-resonant particles (high f2o) in the terminology 
of this paper. This is indeed correct, for the outer-scale 
turbulent fluctuations that we can simulate. If a slow 
wave cascade is present at 1 AU, our /3 = 1 simulations 
in Figure 12 demonstrate that most of the slow- wave en- 
ergy on large-scales is likely lost to particle heating, in 
particular parallel heating of the protons. In addition 
to being important for the thermodynamics of the so- 
lar wind, the slow waves are one of the primary sources 
of density fluctuations in anis otropic Alfvenic turbulence 
(Lithwick & Goldreich 2001); if they are indeed largely 
damped at large scales when /3 ~ 1, this suppresses the 
contribution of slow waves to the small-scale density fluc- 
tuations that are direc tly measured in the solar wind (see 
Chandran et al.||2009 l. 

Although the slow wave component of the cascade can 
lose a signiflcant fraction of its energy on large scales 
to /iV|ji? acceleration of particles (Fig. 12 1, the Alfven- 
wave component of the cascade does not produce signifl- 
cant particle heating until (1) the Alfven wave frequency 
becomes comparable to the cyclotron frequency of parti- 
cles that have a signiflcant density (e.g., protons or he- 
lium) or (2) the perpendicular wavelength of the Alfven 
waves becomes comparable to the proton Larmor radius 
Tl,p, at which point the Alfven w ave cascade transitio ns 
to a kinetic Alfven wave cascade (Howes et al. 2008a). 



From Figure [T0| we conclude that the maximum gy- 
rofrequency of particles that can be cyclotron-resonant 



with the turbulence is fi^ 



2fc 



2/3 



k]iin'"A, where 



\L,max 

we have identifled kmax — k±^max given the anisotropy 
of the Alfvenic cascade. Using the parameters for the 
solar wind above, we estimate that the maximum gy- 
rofrequency of a particle that can be cyclotron-resonant 



Hz <C rip — 0.15 Hz. As a result direct cyclotron res- 
onance is not important at kj__rnax ^ instead, the 
dissipation of the Alfvenic component of weakly com- 
pressible MHD turbulence occurs via the kinetic Alfven 
wave cascade launched when kj_ „iax ^ results 
based on the direct integration of particle orbits in fully 
developed MHD turbulence support previous work that 
reached the same conclusion using linea r theory cascade 
models (e.g., [Qua taert & Gruzinov 1999{ |Cranmer fc van 
[Ballegooijen 20 03f|Howes et al.||2 008a") . 

The dissipation of anisotropic Kinetic Alfven wave tur- 



bulence having fcj^ ^ '^Lp ^^^^^ ^'-'^ fully understood; 
possibilities include Landau resonance (Howes et al. 
2008ai), secondary instabiliti es that gen erate cyclotron 
frequency waves (Markovskii et al. 2006), dissipation in 
current sheets (Drake et al.^2009a) , and stochastic ion or- 
bits created by sutticiently large amplitude waves ( John- 
son & Cheng 120011 iVoitenko & Goossensl 120041). Our 



results, together with the empirical evidence for energet- 
ically important kinetic Alfven w aves in the near-Earth 
solar wind (e.g., Sahraoui et al. 2009) , strongly suggest 
that the anisotropic proton and ion temperatures in the 
solar wind are in part a consequence of heating by this 
kinetic Alfven wave cascade. 
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Fig. 13. — Standard deviation of the magnetic moment /i as a function of the particles' gyrofrequency Qg, after an integration time of 
~ L/cg. After the particle integration begins we either let the turbulence decay (solid line) or keep driving it (dashed line). The parameters 
of the simulation are summarized in Table [l] The continued driving of the turbulence (dashed line) introduces spurious high frequencies 
into the problem which artificially increase me magnetic moment at high f!o • 
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APPENDIX 

A. INTERACTION OF PARTICLES WITH THE TURBULENT DRIVING 



As mentioned in [ 3.2. 1[ driving the turbulence can have undesirable effects on the particles' motion. Our default 
method for driving the turbulence introduces random kicks to the velocity field at each MHD timestep. This artificially 
introduces variations that are much faster than the MHD evolution, and potentially faster than the particles' gyration. 
As a result, the conservation of the magnetic moment can be artificially modified. Figure 13 compares the standard 
deviation of the magnetic moment in two simulations: we either continue driving the turbulence after the particle 
integration is turned on (dashed) or let the turbulence decay (solid). The particles initially have uy q = 1.0 Cs and 
jUo = l-0Cg/i3o, with no variance in ^. After an integration time ~ L/cs, the two simulations produce quite different 
dispersions in /i, in particular at high fioi where the gyrofrequency is much larger than the true frequencies of the 
fiuctuations resolved in the simulation. The continued driving introduces spurious high frequencies into the problem, 
which significantly change the resulting particle heating. This can be remedied by driving with a decorrelation time 
comparable to the outer-scale turnover time of the turbulence, but we choose to be conservative and study the particle 
heating in decaying turbulence. 

B. TESTS 

After writing the particle integrator, we verified it with a series of tests. We ran it with several simple field 
configurations and checked that we recovered known solutions. We also studied how the results depended on the 
chosen time step. Here we summarize some of our tests and their results. 

B.l. Constant magnetic field 

The simplest test is to confirm that one obtains the well-known helical motion in the case of a constant and uniform 
magnetic field and no electric field. For this problem, our method conserves the particles' energy and magnetic moment 
to machine accuracy. The Boris algorithm tends to produce a slightly low gyrofrequency and a slightly large gyroradius. 
However, with our choice of 40 time steps per gyration, the relative error is ~ 10^^. Note also that this simply produces 
a systematic change in the definition of the charge to mass ratio of the particles; there is no secular change in time. 
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Fig. 14. — Standard deviation of the magnetic moment as a function of the particles' gyrofrequency Qq, due to interaction with a 



single parallel-propagating Alfven wave, after an integration time of ~ L/cs. The particles' initial velocities are u 



The curves correspond to different values of /3, but in all cases the wave has the same amplitude {5v 
(fc = 4 X 2tt/L). 
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Fig. 15. — Scatter plot in velocity space {u^^,uj_), for cyclotron-resonant particles {Qg = GOcg/L) interacting with a plane Alfven wave 
{5v = 0.1 Cs; fc = 4 X 2tt/L) at /3 = 0.3. The gray dashed circle has its center at (ha, 0) — with va = 2.6cs - and represents constant energy 
in the frame of the wave. The particles initially all have the same values of and fi; the black dots show the location of the particles after 
interacting with the wave for a time ~ L/cs- As required analytically, the particles move only along the gray dashed curve. 



B.2. Alfven wave 

We let the particles evolve in the presence of a single Alfven wave propagating parallel to the background magnetic 
field. The particles' velocities are ini tiali zed according to the delta function /o from equation ( 17 1, so they all start with 
the same magnetic moment. Figure 14 shows cr(/i) after an integration time of ~ L/cg for ditterent values of f3 (and 
thus different values of va)- The magnetic moment is conserved to high accuracy for particles with gyrofrequencies 
much larger than the frequency of the wave. We also observe several resonances that produce large changes in fi. 
These are consistent with cyclotron resonances, and their positions scale as l3~^^'^, as predicted by linear theory. 

An important feature of Alfven waves is that the perturbed electric field vanishes in the frame of the wave. As a 
consequence, the energy of a particle cannot change in this frame. Back in the frame of the fluid: 



u]_ + (un - VaY 



(Bl) 



where e is the initial energy of the particle. Particles thus evolve on a sphere in velocity space. Figure [15] shows how 
cyclotron-resonant particles are scattered by an Alfven wave in velocity space. The wave has an amplitude dv = 0.1 Cg 
and wavenumber fc = 4 x 2tt/L; (3 — 0.3. The fact that particles remain on a circle centered on va confirms that our 
integrator accurately reproduces the properties of the resonance. 
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B.3. Particle's timestep 

As explained in |3.2.1[ for particles with normalized gyrofrequency flo, our integration time step is chosen to be 

/ 1 27r 1 ^ \ , ^ 

^tparUcle = mm 7^ , -^AtMHD (B2) 

where Ni corresponds to the minimal number of particle time steps per gyration, and N2 corresponds to the minimal 
number of particle time steps within one MHD time step. In our simulations, we use Ni — 40 and N2 = 10. For 
^2^3, there was no significant change in the results. Similarly, so long as A^i > 10, we found that the results were 
converged to better than 1%. 

C. INTERPOLATION METHOD 

This section gives more details on how we interpolate the fields E and B from the MHD grid to the particle's 
position. We use a directionally split interpolation algorithm. The fields are averaged with weights that are the 
product of four 1-dimensional weights. If {x,y,z,t) is the current position of the particle, iq, jo, the indices of the 
nearest cell, and Iq the index of the nearest MHD timestep: 

io + i Jo + 1 feo + 1 io + 1 

^= Y w,{x)wj{y)wk{z)wi{t)F,,j^k,i (CI) 

i=io-l J=Jo-l fc=fco-l /=/o-l 

where the F stands for E or B, and the overline represents interpolation. 

Choosing appropriate 1-D weights is crucial. In particular, ill-chosen expressions may lead to discontinuities in the 
interpolated fields as particles cross from one cell to another, and have proved to produce spurious variation in the 



particles' energy. As explained in ^ 3.2.2 we use the Triangular Shaped Cloud (TSC) method, which ensures that the 
interpolated fields have smoothness m both space and time. In the x direction, for instance: 

where Xi^ is the x position of the nearest cell's center and Ax is the grid spacing. Weights in the j/, z and t directions 
have similar expressions. 



